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pH , Abstract 

o ■ 

Q ■ The surface-impedance matrix method is used to study interfacial 

waves polarized in a plane of symmetry of anisotropic elastic materials. 

Although the corresponding Stroh polynomial is a quartic, it turns out 
^ ' to be analytically solvable in quite a simple manner. A specific appli- 

^vq . cation of the result concerns the calculation of the speed of a Stoneley 

^O , wave, polarized in the common symmetry plane of two rigidly bonded 

anisotropic solids. The corresponding algorithm is robust, easy to im- 
^r . plement, and gives directly the speed (when the wave exists) for any 

orientation of the interface plane, normal to the common symmetry 



plane. Through the examples of the couples (Aluminum)- (Tungsten) 
and (Carbon/epoxy)- (Douglas pine), some general features of a Stone- 
ley wave speed are verified: the wave does not always exist; it is faster 
than the slowest Rayleigh wave associated with the separated half- 
spaces. 

Keywords: Interfacial waves; Bi-material; Surface-impedance matrix; Stroh 
formalism; Anisotropic elasticity. 



1 Introduction 

The semiconductor industry makes great use of wafer bonding, a process 
which allows two different materials to be rigidly and permanently bonded 
along a plane interface, thus producing a composite bi-material [1]. World- 
wide, there are now several hundreds of wafer bonding patents deposited 
yearly [2]. A similar process, fusion bonding, is used by the polymer indus- 
try to bring together two parts of different solid polymers, thus enabling the 
manufacture of a heterogeneous bi-material with specific properties [3, 4]. In 
the first case, wafers of two different crystals are stuck together through van 
der Waals forces, after their surface has been mirror-polished; in the second 
context, a fusion process takes place at the interface, followed by a cooling 
and consolidating period. Whichever the process, it seems important to be 
able to inspect the strength of the bonding, possibly through non-destructive 
ultrasonic evaluation. This is where the study of Stoneley waves (rigid con- 
tact) and of slip waves (sliding contact) is relevant (see Rokhlin et al. [5, 6] 
or Lee and Corbly [7]). In the continuum mechanics literature, great con- 
tributions can be found on the theoretical apprehension of these waves. In 
particular, Barnett et al. [8] for Stoneley waves and Barnett et al. [9] for slip 
waves have provided a rigorous and elegant corpus of results for their possible 
existence and uniqueness, based on the Stroh formalism [10, 11]. However 
very few simple numerical "recipes" exist to compute the speed (and then the 
attenuation factors, partial modes, and profiles) of these waves when they 
exist. 

When the two materials have at least orthorhombic symmetry and their 
crystallographic axes are aligned, the analysis can be conducted in explicit 
form as is best summarized in the article by Chevalier et al. [12]. This explicit 
analysis is possible because for waves proportional to expik{xi + px2 — vt) 
where k, p, v are the respective wave number, attenuation factor, and speed 
of the wave, and xi, X2 (aligned with two common crystallographic axes) 
are the respective directions of propagation and attenuation, the equations 
of motion lead to a propagation condition which is a quadratic in p'^; then 
the relevant roots p can be found in terms of the stiffnesses Cij, C*- for each 
half-space, of the mass densities p, p*, and of the speed v. After construction 
of the general solution to the equations of motion (a linear superposition of 
the partial modes), the boundary condition at the interface (rigid or sliding 
contact) yields the secular equation, of which f is a root. If the crystallo- 
graphic axes of the two orthorhombic materials do not coincide, or when at 




Fig. 1: Plane interface (at X2 = 0) between two different anisotropic materials 
with a common symmetry plane (at Z = 0). 



least one material is monoclinic or triclinic, then the situation becomes much 
more intricate. In the very particular case where the bi-material is made of 
the same material above and below the plane interface, with crystallographic 
axes symmetrically misoriented, a Stoneley wave can be found along the bi- 
sectrix of the misorientation angle [13, 14, 15]. Otherwise, analytical methods 
fail because in general the propagation condition is a sextic in p, unsolvable 
[16] in the Galois sense. 

Now consider the case of a bi-material made of two different anisotropic 
half-spaces with a common symmetry plane Z = 0, orthogonal to the in- 
terface, and with only one common crystallographic axis (the Z axis), see 
Figure 1. Then the in-plane strain decouples from the anti-plane strain [10] 
and in each half-space, the propagation condition factorizes into the prod- 
uct of a quadratic in p (associated with the anti-plane strain) and a quartic 
in p (associated with the in-plane strain). Solving analytically this quartic 
and identifying the two qualifying roots unambiguously in order to be able 
to write down the boundary condition may have seemed a formidable task. 
However, it was recently shown in Fu [17] that such a quartic can in fact be 
solved and the two qualifying roots identified in quite a simple manner, using 
results from algebra. Once the roots p are known, the complete resolution 
of the problem flows out naturally because the surface-impedance matrices 



M{v) and M*{v) for each region are now known explicitly. This knowledge 
gives in turn an explicit secular equation, in the form f[v) = for some func- 
tion /. For Stoneley waves, / is a monotone decreasing function of v, whose 
only zero, when it exists, is the wave speed. This latter property is worth 
emphasizing: therein lies the superiority of the explicit surface-impedance 
matrix method above others based on algebraic manipulations which ulti- 
mately lead to a multitude of secular equations [18, 19] and/or of spurious 
roots [20, 21, 22, 23, 14, 24, 25, 26]. 

This contribution builds on recent advances by Fu and Mielke [27], Mielke 
and Fu [28], and Fu [17], themselves resting on the major works by Barnett 
and Lothe [29], Chadwick and Smith [30], Ting [11], and Mielke and Sprenger 
[31]. We aim at keeping the algebra to a minimum and at delivering a 
numerical recipe giving the wave speed in a robust manner. The reader who 
is keen on implementing such procedures may skip the next two sections to 
jump directly to Section 4, where they are summarized and presented for the 
Rayleigh wave speed and for the Stoneley wave speed. Section 2 recapitulates 
the basic equations of motion and presents the surface-impedance matrices. 
In Section 3, the quartic evoked above is derived, and then explicitly solved 
for the roots which allow for a localization of the wave near the interface. 
Finally, using the "recipes" of Section 4, examples of Rayleigh and Stoneley 
wave speeds computations are presented in Section 5, and the connection is 
made with some numerical results of Chevalier et al. [12]. 

2 Governing equations 

Consider a bi-material made of two distinct anisotropic materials with a 
common symmetry plane, bonded rigidly along a plane interface, 0:2 = say. 
Let p and Cijks be the mass density and elastic stiffnesses of the body below 
{x2 > 0) and p* and Cj*-;,^ be those of the body above (x2 < 0). The Cijks 
and Cj*-fcs are assumed to satisfy the symmetry relations 

Cijks = Cksij = Cjiks, and C^^^^ = Cf^^^ = Cj^j.^, (1) 

and the strong convexity conditions 

Cijksiijiks > 0, C*jj.gC,ijC,ks > 0, V non-zero real symmetric tensors ^. 

(2) 



The strong ellipticity conditions are given by 

CijksViVkljls > 0, Q^ksViVkljls > 0, V non-zero real vectors ry and 7, 

(3) 
and are implied by the strong convexity conditions (2). Let XY Z and X*Y*Z 
be along the crystallographic axes of each material; the X axis (X* axis) 
makes an angle 9 {6*) with the interface, and Z is normal to their common 
symmetry plane and to X2- Finally, let Xi be an axis such that X1X2Z is a 
rectangular coordinate system. See Figure 1. 

In this context, an interfacial wave is a two-component [10] inhomoge- 
neous plane wave, whose propagation is governed by the equations of motion 
for the mechanical displacement u{xi,X2,t) = [ui,U2]'^, 

CijksUk,sj = pUi {X2 > 0), C*ji^^Uk,sj = P*Ui {X2 < 0), (4) 

and which decays away from the interface, 

■u — 7- as X2 ^ ±00. (5) 

Here and henceforward, a comma denotes differentiation with respect to spa- 
tial coordinates and a dot denotes material time derivative. Since the surfaces 
of the upper half-space and of the lower half-space have unit normals (621) 
and (—621), respectively, the traction vectors on these two surfaces are 

U = —Ci2ksUk,s, ^^ = C*2ksUk,s, "^ = 1,2, (6) 

and by (5), they also decay, 

tj — > as X2 — > +00, t* ^ as X2 — )■ —00. (7) 

Without loss of generality, the interfacial wave is assumed to propagate along 
the xi-direction and to have unit wave number, so that 

u = z{iX2)e'^'''-'''^ + c.c., (8) 

where v is the propagation speed and "c.c." denotes the complex conjugate 
of the preceding term. 

Substituting (8) into (4) and (6)1 gives 

Tz" + {R + R^)z' + (g - pv'^I)z = 0, X2> 0, (9) 



and 
where 



-il{ix2)e 



i(xi~vt) 



+ C.C, 



l = Tz' + R^z, 



(10) 



:ii) 



a prime signifies differentiation with respect to the argument ix2, and the 
2x2 matrices T, R, Q are defined by their components, 



T,, 



ik 



a 



i2k2, 



Rik — L/,: 



ilfc25 



Q 



ik 



a. 



ilkl- 



(12) 



Identical results apply for the upper half-space X2 < Q, with each quantity 
replaced by its starred counterpart, except that t* is given by 



t* = ir{ix2)e 



i(xi—vt) 



+ C.C. 



(13) 



Note that satisfaction of the strong ellipticity conditions (3) ensures that T, 
T*, and Q, Q* are all positive definite and hence invertible. 

The surface-impedance matrices M{v) and M*{v) are defined by 



-zZ(0) = M{v)z{0), zZ*(0) = M*{v)z*{0). 



(14) 



In the Stroh [10] formulation, the second-order differential equation (9) is 
written as a system of first-order differential equations for the variables z 
and I. Thus, for the lower half-space (and similarly for the upper half-space). 



^' = iV^, where $, 



N 



Ni N2 

Ns + pv^I Nf 



and 



iVi 



-T-^R^, 



N2 = T-\ iVg 



RT-^R^ 



Q. 



(15) 



(16) 



Here / is the 2x2 identity matrix. 

An important property of the surface impedance matrix is that it is in- 
dependent of the depth, that is 



-l{ix2) = Mz{ix2), 



(17) 



see Ingebrigsten and Tonning [32]. On substituting (17) into (15) and elimi- 
nating z', we obtain 



{(M - iR)T-\M + iR^) -Q + pv^} z{ix2) = 0. 
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(18) 



Since z{ix2) is arbitrary, it then follows that 

(M - iR)T^\M + iR^) -Q- pv^I = 0. 



(19) 



This simple matrix equation satisfied by M was seemingly first derived by 
Biryukov [33], and later rediscovered by Fu and Mielke [27] where it was 
shown how this equation could be exploited to compute the surface speed. 

The general solution of (15) can be constructed by first looking for a 
partial mode solution of the form 



^{1x2] 



^ipX2 



c, 



(20) 



say, where p is a constant scalar (attenuation factor) and <^ is a constant 
vector to be determined. Substituting (20) into (15) leads to the eigenvalue 
problem (A^ — p/4x4)C = 0. Under the assumption of strong ellipticity, the 
eigenvalues of A^ appear as two pairs of complex conjugates when v = 
and they remain so until v = Vc, where Vc is referred to as the limiting 
speed. In this paper we are only concerned with the subsonic case, for which 
< V < Vc- The solution (20) decays as a;2 — )■ 00 only if the imaginary 
part of p is positive. Thus, denoting by pi, p2 the two eigenvalues of N with 
positive imaginary parts, and by C , C the corresponding eigenvectors, a 
general decaying solution is given by 



^{1x2] 



cie 



ipiX2/-il) 



C'> + C2e 



ip2X2/-i2) 



(21) 



where ci, C2 are disposable constants. Hence at the interface. 



^(0) = c,C^ + C2C 



.(2) 



(22) 



where the 2x2 matrices A, B and the column vector c are defined by 



^(l)|^{2) 



It follows from (22) that 

-«/(0) = -iBc = -tBA-^z{0), 

and so from (14) that 

M{v) = -iBA-~\ 



(23) 

(24) 
(25) 
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which is a well-known representation of the surface-impedance matrix. As 
is recalled later in the paper, the surface-impedance matrices are crucial to 
the determination of the interfacial wave speed. In fact, if their explicit 
expressions are found, then an exact secular equation (an equation of which 
the wave speed is the only zero) is also found explicitly. 

3 Explicit expressions for the surface-impe- 
dance matrices 

Here it is seen that the relevant roots to the characteristic equation det(A^ — 
plixi) = can be obtained explicitly, without the uncertainty which one 
encounters in, for instance, solving a cubic for p^. Indeed, because the char- 
acteristic equation is a quartic in p, the qualifying roots are those with a 
positive (negative) imaginary part to ensure decay with distance from the 
interface in the lower (upper) half-space. Thus, for the lower half-space say, 
they must be of the form ai + iPi, 0:2 + iP2, where (3i and /32 are positive 
real numbers. This situation is in sharp contrast with the case of a wave 
propagating in a symmetry plane; then the characteristic equation is a cu- 
bic in p^ and the relevant roots for the lower half-space can come in one of 
two forms: either as iai, ia2, ia^, where ai, 0:2, 0:3 are positive, or as iai, 
±a + if3, where ai and f3 are positive. Although the roots of a cubic are 
seemingly easier to obtain analytically than those of a quartic, determining 
which of these two forms applies is a tricky matter. Here, once pi, p2 and pi, 
P2 are known, M{v) and M*{v) can be constructed and the interfacial wave 
speed can be computed directly. The analysis below is conducted for the 
lower half-space {x2 > 0) and indications are given at the end of the section 
on how to adapt it to the upper half-space. 

First write the four- component vectors Q^' and Q"^' as 



^(1) 






c 



(2) 



5(2) 



(26) 



The vectors a^^\ a'^'^^ are determined from (9), that is from 

\plT + pk{R + R^) + Q-pv^I]a^''^ =0, A; = 1,2, (27) 

and the vectors 6*- •*, b^ ', are computed from (11) according to 

b« = (pfcT + R^)a^^\ k = l,2. (28) 



It then follows from (23) that 
where 



fed) 1^(2) 



Pi 

V2 



T 



TAQ + W A, 



(29) 



(30) 



The two eigenvalues pi and p2 are determined from det (A^ — pi) = 0, or 
equivalently, from 



det [/T + p{R + R^) + Q- pv^I] = 0. 



(31) 



This characteristic equation, called the propagation condition, is a quartic 
which can be written as 



p^ + dsp^ + d2P^ + dip + (io = 0, 



(32) 



say, where dQ,di,d2, d^ are real constants. 

The usual strategy for solving a quartic equation such as (32) is to use 
a substitution to eliminate the cubic term; see Bronshtein and Semendyayev 
[34]. Thus, with the substitution p = q — d^/i, equation (32) reduces to 

q'^ + rq'^ + sq + h = 0, (33) 

where 

r = d2-ldl, s = (ii- 1^2(^3 + 1^3, h = do~\dids + -^d2dl~^dl. (34) 

The behavior of the roots of (33) depends on the cubic resolvent 



z^ + 2rz^ + (r^ - Ah)z - s" 



0. 



(35) 



In particular, Eq.(33), and hence Eq.(32), have two pairs of complex con- 
jugate solutions if and only if (35) has three real roots zi, Z2, z^ such that 
-22 < -23 < < 2;i. Then, 



d^ 

Pi = |[sign(s)v/ir+ «(v/^i2 + v^^is)] - -j, 

P2 = l[-sign{s)y^+^{^/^z^-^/^z^)]- ^, 



d:^ 



(36) 



where sign(s) equals 1 if s is non-negative and —1 otherwise (this definition 
overrides the standard definition in which sign(O) = 0). The three roots of 
(35) are 

2A^/^cos(0)-|r, 2Ai/3cos(0 + 27r/3)-|r, 2Ai/3cos(0+47r/3)-|r, (37) 

where 

X = ^{12h + ry/\ cos30 = f (12/i + r2)-3/2(^r3 + s2_|^/,). (33) 

Without loss of generality it can be assumed that < 30 < vr. It is easy to 
show that in this interval, 

cos(0 + 27r/3) < cos(0 + An/S) < cos 0. (39) 

It then follows that the three roots of (35) can explicitly be identified as 

zi =2Ai/3cos(0)-|r, 

Z2 = 2A^/^ cos(0 + 27r/3) - |r, 

Z3 = 2A1/3 cos(0 + 47r/3) - |r. (40) 

Since Z2 7^ z^, we may further deduce that 7^ so that < < n/3. 

Now all the ingredients are in place to compute explicitly pi, p2', then A, 
B; and ultimately, M{v), which is all that is required to find the speed of the 
Rayleigh wave propagating at the interface between the lower half-space and 
a vacuum. To compute the speed of Stoneley waves and slip waves, M* (v) is 
needed. For the upper half-space X2 < 0, the analysis above can be repeated 
by starring all the quantities involved. Now the qualifying roots pi, P2 must 
have negative imaginary parts in order to satisfy the decaying condition and 
so (36) is replaced by 

pl = |[sign(s*)v/if - z(v/^+ v^^if)] - -f, 

pl = |[-sign(s*)v/ir- ^(v/^- 7=4)] - f • (41) 

As a result, M*{v) can be obtained from M{v) by first replacing the ma- 
terial constants by their starred counterparts and then taking the complex 
conjugate. 
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4 How to find the wave speed 

This section presents simple algorithms that can be used to compute directly 
the surface wave speed and the interfacial wave speed. They rely on the use 
of a symbolic manipulation package such as Mathematica. The algorithm 
for Rayleigh waves (solid/vacuum interface) is presented in detail, and then 
modified for Stoneley waves (rigid solid/solid interface). 

4.1 Rayleigh waves 

A Rayleigh wave satisfies the traction free boundary condition Z(0) = 0. 
Thus, from (14)i, its speed is given by 

det M{v) = det {-zBA^^) = -det {TAnA'^ + R^) = 0. (42) 

Given a set of material constants CijM and p, the unique surface-wave 
speed is found by using the following robust numerical procedure: 

(i) Enter the values of the stiffnesses into the definitions (12) of T, R, Q, 
that is 



Q 



Ciiii Cni2 



R 



C1II2 C\Y2,2 
C\2\2 C22\2 



T 



C1212 (^1222 

C1222 C*2222 



m 



IV 



fvi 



Expand the quartic (31) in p, and obtain the coefficients d\^ ^2, (is by 
comparing it to (32). 

Enter the coefficients r, s, h according to (34) and then enter pi,p2 
according to (36) and (40). 

Define a'^^^ and a*^^-* according to 



a' 



PiTi2+Pfc(i?l2 + i?2l)+Q2l' 

-vlTxx - '2,pkRii - Qii + pv^ 



k = 1,2. 



(43) 



Define A from (29) 1, Q from (30). 



Use (42) and the command FindRoot in Mathematica to solve det M{v) 
0. 



To facilitate calculations, the term sign{s) ^/z^ in (36) may be replaced 
by s/y/z2Z3. If necessary, the solution for M can be checked by substituting 
the result into the matrix equation (19). 
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4.2 Stoneley waves 

A Stoneley wave must satisfy the continuity conditions 

z{0) = z*{0), l{0) = l*{0). 

It then follows from (14) that [M{v) + M*{v)]z{0) = and so the Stoneley 
wave speed is determined by 

det[M{v) + M*{v)] = 0, (44) 

where 

M{v) = -i[TAQA-^ + R^], M*{v) = i[T*A*Q*{A*)-'^ + {R*f]. (45) 

Given two sets of material constants Cijki, C*j^i and p, p*, the unique 
Stoneley wave speed, if it exists, is found by the following robust numerical 
procedure: 

• Follow steps (i) to (iv) of the algorithm described in the preceding 
subsection twice: once for the lower half-space, and once for the upper 
half-space, replacing each quantity but v by their starred counterpart, 
taking care in step (iii) that p\, p^ are defined by (41). 

• Use (44), (45) and the command FindRoot in Mathematica to find the 
speed of the Stoneley wave, when it exists. If Mathematica is unable 
to find a root, then the Stoneley wave does not exist. 

5 Examples 

We now apply the algorithms to specific materials. We use the data and 
results of Chevalier et al. [12], where ^ = ^* = 0, as a guideline. By varying 
these angles, we find that there are situations where a Stoneley wave does not 
exist and that, when it does exist, it is faster than the slower Rayleigh wave 
associated with either of the two half-spaces; these features were proved for 
any anisotropic crystal by Barnett et al. [8]. We also compute the smallest 
imaginary part of the attenuation coefficients, which is Q{p2) according to 
(36); this quantity is related to the penetration depth of the interfacial wave 
into the substrates: the smaller it is, the deeper is the penetration. 
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Fig. 2: Rayleigh (thin curves) and Stoneley (thick curve) wave speeds for 
a bi-material made of aluminum (rotated Y-cnt about the Z axis) and of 
tungsten (symmetrically rotated Y-cut about the Z axis). 



5.1 (Aluminum)- (Tungsten) bi-material 

For a bi-material made of aluminum above {x2 < 0) and of tungsten below 
{x2 > 0), Chevalier et al. [12] found a Stoneley wave propagating at speed 
2787 m/s when 6 = 6* = 0. We recovered this result and extended it to the 
consideration of the bi-material obtained when the half-spaces are rotated 
symmetrically about the Z axis before they are cut and bonded, 



0° < ^ < 90°, 9* 



-9. 



(46) 



Figure 2 shows that the Stoneley wave (thick curve) exists for all angles; it 
is faster than the Rayleigh wave for a half-space made of tungsten (lower thin 
curve) and slower than the Rayleigh wave for a half-space made of aluminum 
(upper thin curve). 

Figure 3 displays the smallest imaginary part of the attenuation factors for 
each wave. We find that for the Stoneley wave, this quantity is intermediate 
between the corresponding quantities for the Rayleigh waves, indicating a 
similar localization. 
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Fig. 3: Attenuation factors for Rayleigh (thin curves) and Stoneley (thick 
curve) wave speeds for a bi-material made of aluminum (rotated F-cut about 
the Z axis) and of tungsten (symmetrically rotated F-cut about the Z axis). 



5.2 (Douglas pine)-(Carbon/epoxy) bi-material 

For a bi-material made of carbon/epoxy above (x2 < 0) and of Douglas pine 
below (x2 > 0), Chevalier et al. [12] found a Stoneley wave propagating 
at speed 1353.7 m/s when 6* = and 6 = 90°. Here we investigate what 
happens to this wave when the half-space below is rotated about the Z axis 
before it is cut and bonded, while the half-space above is left untouched, 



0° < ^ < 90°, 9* = 0. 



(47) 



Figure 4 shows that the Stoneley wave indeed exists at ^ = 90° and 
in the neighborhood of that angle, approximatively in the range: 72.4° < 
6 < 107.6°. In that range, the Rayleigh wave for a half-space made of car- 
bon/epoxy cut at an angle 6 (thin varying curve) is slower than the Rayleigh 
wave for a half-space made of Douglas pine cut at an angle 6* = (thin 
horizontal curve). The Stoneley wave, when it exists (thick curve), is al- 
ways faster than the former, and either faster (for 72.4° < 9 < 75.3° and for 
104.7° <e < 107.6°) or slower (for 75.3° <e < 104.7°) than the latter. 

Finally, Figure 5 displays the smallest imaginary part of the attenuation 
factors for each wave. We find that this quantity is for the Stoneley wave 
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Fig. 4: Rayleigh (thin curves) and Stoneley (thick curve) wave speeds (km/s) 
for a bi-material made of Douglas pine (F-cut) and of carbon/epoxy (rotated 
Y-cut about the Z axis). 



Q{q) 



0.02 



Fig. 5: Attenuation factors for Rayleigh (thin curves) and Stoneley (thick 
curve) wave speeds for a bi-material made of Douglas pine (F-cut) and of 
carbon/epoxy (rotated Y-cut about the Z axis), in the range of common 
existence. 

between one sixth and one tenth of that for the Rayleigh waves, indicating 
a much deeper penetration. 
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